{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "14\n"
     ]
    }
   ],
   "source": [
    "import sys,os\n",
    "add_path = os.path.abspath('../python')\n",
    "if not add_path in sys.path:\n",
    "    sys.path.append(add_path)\n",
    "print(sys.path.index(add_path))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import time\n",
    "import numpy as np\n",
    "from numpy import vstack, hstack\n",
    "from segbench import csaAssign, correspondPixels,correspondVoxels, matchEdgeMaps2D, matchEdgeMaps3D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 1.  5.  1.]\n",
      " [ 2.  6.  1.]\n",
      " [ 3.  4.  5.]]\n"
     ]
    }
   ],
   "source": [
    "data = np.array([[1, 4, 3],\n",
    "    [2, 5, 3],\n",
    "    [3, 6, 3],\n",
    "    [1, 5, 1],\n",
    "    [2, 6, 1],\n",
    "    [3, 4, 5]])\n",
    "n = 6\n",
    "e1 = csaAssign(n,data)\n",
    "print e1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "6.949e-02s\n",
      "1000 500130 [  4.00000000e+00   1.99300000e+03   1.00000000e+00]\n"
     ]
    }
   ],
   "source": [
    "# Unit test for csaAssign.\n",
    "n = 1000\n",
    "dg = np.random.rand(n,n)>.5\n",
    "m = np.sum(dg.flat)\n",
    "i = np.array(dg.ravel().nonzero()).astype(float)\n",
    "g1 = 1+np.floor(np.divide(i,n))\n",
    "g2 = (1+np.floor(np.divide(i,n)),\n",
    "     1+ np.mod(i,n)+n,\n",
    "     1+np.floor(np.random.rand(1,m))*1000)\n",
    "g2 = np.array(np.vstack(g2)).T\n",
    "t = time.time()\n",
    "e2 = csaAssign(n*2,g2)\n",
    "print '{:.3e}s'.format(time.time() - t)\n",
    "esum = np.sum(e2,axis=0)\n",
    "if not all(esum[0:2]==(n*(n+1)/2,n*(n+1)/2+n*n)): raise ValueError(\"bug\")\n",
    "if not np.sum(esum[0:2])==(2*n*(2*n+1)/2): raise ValueError(\"bug\")\n",
    "print n,m, e2[3][:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "n=1 thresh=0.00e+00 1.32e-03s 2 0 2\n",
      "n=1 thresh=1.11e-01 5.70e-05s 2 0 2\n",
      "n=1 thresh=2.22e-01 3.10e-05s 1 0 1\n",
      "n=1 thresh=3.33e-01 5.20e-05s 2 0 2\n",
      "n=1 thresh=4.44e-01 3.00e-05s 1 0 1\n",
      "n=1 thresh=5.56e-01 3.10e-05s 1 0 1\n",
      "n=1 thresh=6.67e-01 2.91e-05s 0 0 0\n",
      "n=1 thresh=7.78e-01 2.81e-05s 1 0 1\n",
      "n=1 thresh=8.89e-01 3.50e-05s 0 0 0\n",
      "n=1 thresh=1.00e+00 4.89e-05s 0 0 0\n",
      "n=2 thresh=0.00e+00 9.80e-05s 8 6 2\n",
      "n=2 thresh=1.11e-01 1.81e-04s 8 6 2\n",
      "n=2 thresh=2.22e-01 1.53e-04s 8 6 2\n",
      "n=2 thresh=3.33e-01 5.98e-05s 5 2 3\n",
      "n=2 thresh=4.44e-01 5.91e-05s 5 2 3\n",
      "n=2 thresh=5.56e-01 4.91e-05s 4 2 2\n",
      "n=2 thresh=6.67e-01 3.10e-05s 2 0 2\n",
      "n=2 thresh=7.78e-01 3.10e-05s 3 0 3\n",
      "n=2 thresh=8.89e-01 5.58e-05s 0 0 0\n",
      "n=2 thresh=1.00e+00 3.29e-05s 0 0 0\n",
      "n=3 thresh=0.00e+00 2.11e-04s 18 16 2\n",
      "n=3 thresh=1.11e-01 2.20e-04s 16 12 4\n",
      "n=3 thresh=2.22e-01 2.09e-04s 13 8 5\n",
      "n=3 thresh=3.33e-01 2.47e-04s 14 12 2\n",
      "n=3 thresh=4.44e-01 9.82e-05s 7 2 5\n",
      "n=3 thresh=5.56e-01 8.80e-05s 7 2 5\n",
      "n=3 thresh=6.67e-01 1.01e-04s 5 2 3\n",
      "n=3 thresh=7.78e-01 8.80e-05s 5 2 3\n",
      "n=3 thresh=8.89e-01 3.19e-05s 2 0 2\n",
      "n=3 thresh=1.00e+00 3.10e-05s 0 0 0\n",
      "n=5 thresh=0.00e+00 5.65e-04s 50 48 2\n",
      "n=5 thresh=1.11e-01 4.20e-04s 44 36 8\n",
      "n=5 thresh=2.22e-01 3.05e-04s 41 30 11\n",
      "n=5 thresh=3.33e-01 1.97e-04s 33 16 17\n",
      "n=5 thresh=4.44e-01 1.27e-04s 27 10 17\n",
      "n=5 thresh=5.56e-01 9.11e-05s 22 6 16\n",
      "n=5 thresh=6.67e-01 6.20e-05s 15 4 11\n",
      "n=5 thresh=7.78e-01 5.51e-05s 13 2 11\n",
      "n=5 thresh=8.89e-01 6.10e-05s 8 4 4\n",
      "n=5 thresh=1.00e+00 8.08e-05s 0 0 0\n",
      "n=8 thresh=0.00e+00 1.36e-03s 128 126 2\n",
      "n=8 thresh=1.11e-01 1.04e-03s 113 100 13\n",
      "n=8 thresh=2.22e-01 7.99e-04s 97 76 21\n",
      "n=8 thresh=3.33e-01 7.49e-04s 87 58 29\n",
      "n=8 thresh=4.44e-01 3.67e-04s 74 36 38\n",
      "n=8 thresh=5.56e-01 1.69e-04s 50 16 34\n",
      "n=8 thresh=6.67e-01 1.48e-04s 50 12 38\n",
      "n=8 thresh=7.78e-01 8.42e-05s 38 6 32\n",
      "n=8 thresh=8.89e-01 7.20e-05s 17 4 13\n",
      "n=8 thresh=1.00e+00 4.41e-05s 0 0 0\n",
      "n=13 thresh=0.00e+00 5.27e-03s 338 332 6\n",
      "n=13 thresh=1.11e-01 3.66e-03s 311 282 29\n",
      "n=13 thresh=2.22e-01 3.05e-03s 261 212 49\n",
      "n=13 thresh=3.33e-01 2.22e-03s 229 152 77\n",
      "n=13 thresh=4.44e-01 1.71e-03s 201 108 93\n",
      "n=13 thresh=5.56e-01 1.16e-03s 143 62 81\n",
      "n=13 thresh=6.67e-01 5.94e-04s 105 36 69\n",
      "n=13 thresh=7.78e-01 2.32e-04s 68 10 58\n",
      "n=13 thresh=8.89e-01 2.14e-04s 46 12 34\n",
      "n=13 thresh=1.00e+00 8.61e-05s 0 0 0\n",
      "n=22 thresh=0.00e+00 1.81e-02s 968 966 2\n",
      "n=22 thresh=1.11e-01 1.35e-02s 867 774 93\n",
      "n=22 thresh=2.22e-01 8.52e-03s 737 554 183\n",
      "n=22 thresh=3.33e-01 7.56e-03s 652 446 206\n",
      "n=22 thresh=4.44e-01 3.29e-03s 520 262 258\n",
      "n=22 thresh=5.56e-01 2.58e-03s 419 174 245\n",
      "n=22 thresh=6.67e-01 1.28e-03s 314 90 224\n",
      "n=22 thresh=7.78e-01 5.92e-04s 201 44 157\n",
      "n=22 thresh=8.89e-01 2.89e-04s 125 14 111\n",
      "n=22 thresh=1.00e+00 2.32e-04s 0 0 0\n",
      "n=36 thresh=0.00e+00 5.51e-02s 2592 2584 8\n",
      "n=36 thresh=1.11e-01 4.23e-02s 2306 2052 254\n",
      "n=36 thresh=2.22e-01 3.22e-02s 2008 1556 452\n",
      "n=36 thresh=3.33e-01 2.20e-02s 1706 1114 592\n",
      "n=36 thresh=4.44e-01 1.55e-02s 1419 756 663\n",
      "n=36 thresh=5.56e-01 8.03e-03s 1126 486 640\n",
      "n=36 thresh=6.67e-01 3.58e-03s 835 260 575\n",
      "n=36 thresh=7.78e-01 1.47e-03s 569 114 455\n",
      "n=36 thresh=8.89e-01 4.72e-04s 271 26 245\n",
      "n=36 thresh=1.00e+00 1.96e-04s 0 0 0\n",
      "n=60 thresh=0.00e+00 2.03e-01s 7200 7184 16\n",
      "n=60 thresh=1.11e-01 1.48e-01s 6436 5760 676\n",
      "n=60 thresh=2.22e-01 1.15e-01s 5579 4290 1289\n",
      "n=60 thresh=3.33e-01 7.26e-02s 4778 3112 1666\n",
      "n=60 thresh=4.44e-01 5.63e-02s 4025 2242 1783\n",
      "n=60 thresh=5.56e-01 3.44e-02s 3231 1426 1805\n",
      "n=60 thresh=6.67e-01 1.41e-02s 2399 776 1623\n",
      "n=60 thresh=7.78e-01 6.06e-03s 1596 348 1248\n",
      "n=60 thresh=8.89e-01 1.63e-03s 825 76 749\n",
      "n=60 thresh=1.00e+00 4.84e-04s 0 0 0\n",
      "n=100 thresh=0.00e+00 9.22e-01s 20000 19942 58\n",
      "n=100 thresh=1.11e-01 4.89e-01s 17741 15682 2059\n",
      "n=100 thresh=2.22e-01 3.80e-01s 15536 12034 3502\n",
      "n=100 thresh=3.33e-01 2.49e-01s 13221 8760 4461\n",
      "n=100 thresh=4.44e-01 1.55e-01s 11036 6158 4878\n",
      "n=100 thresh=5.56e-01 9.49e-02s 8771 3810 4961\n",
      "n=100 thresh=6.67e-01 4.85e-02s 6572 2122 4450\n",
      "n=100 thresh=7.78e-01 1.98e-02s 4491 1022 3469\n",
      "n=100 thresh=8.89e-01 4.28e-03s 2297 248 2049\n",
      "n=100 thresh=1.00e+00 1.10e-03s 0 0 0\n"
     ]
    }
   ],
   "source": [
    "# Unit test for correspondPixels.\n",
    "for n in np.round(np.logspace(0,2,10)).astype(int):\n",
    "    for thresh in np.linspace(0,1,10):\n",
    "        print 'n={:d} thresh={:1.2e}'.format(n,thresh),\n",
    "        bmap1 = (np.random.rand(n,n)>thresh).astype(float)\n",
    "        bmap2 = (np.random.rand(n,n)>thresh).astype(float)\n",
    "        t = time.time()\n",
    "        cost,oc,match1,match2 = correspondPixels(bmap1,bmap2,.005)\n",
    "        print '{:1.2e}s'.format(time.time() - t),\n",
    "        npix = len(bmap1.nonzero()[0]) + len(bmap2.nonzero()[0])\n",
    "        ncor = len(match1.nonzero()[0]) + len(match1.nonzero()[0])\n",
    "        nout = len(((bmap1>0)&(match1==0)).nonzero()[0]) + len(((bmap2>0)&(match2==0)).nonzero()[0])\n",
    "        print npix,ncor,nout\n",
    "        if ncor+nout!=npix: raise ValueError(\"bug\")\n",
    "\n",
    "        match1f = match1.ravel()\n",
    "        idx1 = match1f.astype(int)\n",
    "        a1 = match1f.nonzero()[0]\n",
    "        a2 = match1f[idx1.nonzero()[0]].ravel()\n",
    "        a = vstack((a1,a2)).T\n",
    "\n",
    "        match2f = match2.ravel()\n",
    "        idx2 = match2f.astype(int)\n",
    "        b1 = match2f.nonzero()[0]\n",
    "        b2 = match2f[idx2.nonzero()[0]].ravel()\n",
    "        b = vstack((b2,b1)).T\n",
    "        b = b[b[:,0].argsort(),]\n",
    "        if np.sum(a-b): raise ValueError(\"bug\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "n=1 thresh=0.00e+00 1.02e-04s 0 + 2 = 2 == 2\n",
      "n=1 thresh=1.11e-01 5.60e-05s 0 + 2 = 2 == 2\n",
      "n=1 thresh=2.22e-01 2.88e-05s 0 + 1 = 1 == 1\n",
      "n=1 thresh=3.33e-01 5.29e-05s 0 + 2 = 2 == 2\n",
      "n=1 thresh=4.44e-01 3.70e-05s 0 + 1 = 1 == 1\n",
      "n=1 thresh=5.56e-01 5.70e-05s 0 + 2 = 2 == 2\n",
      "n=1 thresh=6.67e-01 2.50e-05s 0 + 0 = 0 == 0\n",
      "n=1 thresh=7.78e-01 2.91e-05s 0 + 0 = 0 == 0\n",
      "n=1 thresh=8.89e-01 2.79e-05s 0 + 0 = 0 == 0\n",
      "n=1 thresh=1.00e+00 3.10e-05s 0 + 0 = 0 == 0\n",
      "n=2 thresh=0.00e+00 9.20e-05s 6 + 2 = 8 == 8\n",
      "n=2 thresh=1.11e-01 6.60e-05s 6 + 1 = 7 == 7\n",
      "n=2 thresh=2.22e-01 7.01e-05s 2 + 4 = 6 == 6\n",
      "n=2 thresh=3.33e-01 6.01e-05s 4 + 1 = 5 == 5\n",
      "n=2 thresh=4.44e-01 5.39e-05s 2 + 2 = 4 == 4\n",
      "n=2 thresh=5.56e-01 3.31e-05s 0 + 3 = 3 == 3\n",
      "n=2 thresh=6.67e-01 2.98e-05s 0 + 2 = 2 == 2\n",
      "n=2 thresh=7.78e-01 5.48e-05s 4 + 1 = 5 == 5\n",
      "n=2 thresh=8.89e-01 3.10e-05s 0 + 2 = 2 == 2\n",
      "n=2 thresh=1.00e+00 2.60e-05s 0 + 0 = 0 == 0\n",
      "n=3 thresh=0.00e+00 2.22e-04s 16 + 2 = 18 == 18\n",
      "n=3 thresh=1.11e-01 1.72e-04s 12 + 4 = 16 == 16\n",
      "n=3 thresh=2.22e-01 1.40e-04s 8 + 4 = 12 == 12\n",
      "n=3 thresh=3.33e-01 5.91e-05s 2 + 7 = 9 == 9\n",
      "n=3 thresh=4.44e-01 9.20e-05s 6 + 6 = 12 == 12\n",
      "n=3 thresh=5.56e-01 6.41e-05s 4 + 5 = 9 == 9\n",
      "n=3 thresh=6.67e-01 5.70e-05s 2 + 5 = 7 == 7\n",
      "n=3 thresh=7.78e-01 3.79e-05s 0 + 5 = 5 == 5\n",
      "n=3 thresh=8.89e-01 3.50e-05s 0 + 3 = 3 == 3\n",
      "n=3 thresh=1.00e+00 2.91e-05s 0 + 0 = 0 == 0\n",
      "n=5 thresh=0.00e+00 5.15e-04s 48 + 2 = 50 == 50\n",
      "n=5 thresh=1.11e-01 3.32e-04s 32 + 10 = 42 == 42\n",
      "n=5 thresh=2.22e-01 3.32e-04s 32 + 9 = 41 == 41\n",
      "n=5 thresh=3.33e-01 1.58e-04s 12 + 17 = 29 == 29\n",
      "n=5 thresh=4.44e-01 2.32e-04s 18 + 12 = 30 == 30\n",
      "n=5 thresh=5.56e-01 7.51e-05s 6 + 16 = 22 == 22\n",
      "n=5 thresh=6.67e-01 6.89e-05s 4 + 12 = 16 == 16\n",
      "n=5 thresh=7.78e-01 6.79e-05s 4 + 8 = 12 == 12\n",
      "n=5 thresh=8.89e-01 3.81e-05s 0 + 5 = 5 == 5\n",
      "n=5 thresh=1.00e+00 3.70e-05s 0 + 0 = 0 == 0\n",
      "n=8 thresh=0.00e+00 1.37e-03s 126 + 2 = 128 == 128\n",
      "n=8 thresh=1.11e-01 1.18e-03s 92 + 16 = 108 == 108\n",
      "n=8 thresh=2.22e-01 8.97e-04s 80 + 23 = 103 == 103\n",
      "n=8 thresh=3.33e-01 7.49e-04s 74 + 22 = 96 == 96\n",
      "n=8 thresh=4.44e-01 3.88e-04s 30 + 37 = 67 == 67\n",
      "n=8 thresh=5.56e-01 2.88e-04s 26 + 31 = 57 == 57\n",
      "n=8 thresh=6.67e-01 1.83e-04s 12 + 27 = 39 == 39\n",
      "n=8 thresh=7.78e-01 7.61e-05s 2 + 25 = 27 == 27\n",
      "n=8 thresh=8.89e-01 4.60e-05s 0 + 5 = 5 == 5\n",
      "n=8 thresh=1.00e+00 4.41e-05s 0 + 0 = 0 == 0\n",
      "n=13 thresh=0.00e+00 4.59e-03s 336 + 2 = 338 == 338\n",
      "n=13 thresh=1.11e-01 3.39e-03s 258 + 37 = 295 == 295\n",
      "n=13 thresh=2.22e-01 2.09e-03s 194 + 61 = 255 == 255\n",
      "n=13 thresh=3.33e-01 1.46e-03s 140 + 71 = 211 == 211\n",
      "n=13 thresh=4.44e-01 1.31e-03s 104 + 95 = 199 == 199\n",
      "n=13 thresh=5.56e-01 7.67e-04s 64 + 77 = 141 == 141\n",
      "n=13 thresh=6.67e-01 4.49e-04s 50 + 77 = 127 == 127\n",
      "n=13 thresh=7.78e-01 1.55e-04s 10 + 60 = 70 == 70\n",
      "n=13 thresh=8.89e-01 8.70e-05s 4 + 23 = 27 == 27\n",
      "n=13 thresh=1.00e+00 5.89e-05s 0 + 0 = 0 == 0\n",
      "n=22 thresh=0.00e+00 1.54e-02s 966 + 2 = 968 == 968\n",
      "n=22 thresh=1.11e-01 1.43e-02s 770 + 93 = 863 == 863\n",
      "n=22 thresh=2.22e-01 9.89e-03s 570 + 180 = 750 == 750\n",
      "n=22 thresh=3.33e-01 6.16e-03s 412 + 220 = 632 == 632\n",
      "n=22 thresh=4.44e-01 5.09e-03s 290 + 243 = 533 == 533\n",
      "n=22 thresh=5.56e-01 2.37e-03s 186 + 228 = 414 == 414\n",
      "n=22 thresh=6.67e-01 1.19e-03s 98 + 227 = 325 == 325\n",
      "n=22 thresh=7.78e-01 5.31e-04s 36 + 177 = 213 == 213\n",
      "n=22 thresh=8.89e-01 3.36e-04s 20 + 86 = 106 == 106\n",
      "n=22 thresh=1.00e+00 1.10e-04s 0 + 0 = 0 == 0\n",
      "n=36 thresh=0.00e+00 5.96e-02s 2584 + 8 = 2592 == 2592\n",
      "n=36 thresh=1.11e-01 4.80e-02s 2070 + 248 = 2318 == 2318\n",
      "n=36 thresh=2.22e-01 2.95e-02s 1556 + 448 = 2004 == 2004\n",
      "n=36 thresh=3.33e-01 2.09e-02s 1118 + 603 = 1721 == 1721\n",
      "n=36 thresh=4.44e-01 1.27e-02s 804 + 638 = 1442 == 1442\n",
      "n=36 thresh=5.56e-01 8.14e-03s 512 + 662 = 1174 == 1174\n",
      "n=36 thresh=6.67e-01 4.29e-03s 292 + 583 = 875 == 875\n",
      "n=36 thresh=7.78e-01 1.90e-03s 144 + 489 = 633 == 633\n",
      "n=36 thresh=8.89e-01 5.52e-04s 32 + 274 = 306 == 306\n",
      "n=36 thresh=1.00e+00 2.17e-04s 0 + 0 = 0 == 0\n",
      "n=60 thresh=0.00e+00 1.89e-01s 7182 + 18 = 7200 == 7200\n",
      "n=60 thresh=1.11e-01 1.31e-01s 5632 + 760 = 6392 == 6392\n",
      "n=60 thresh=2.22e-01 1.13e-01s 4354 + 1265 = 5619 == 5619\n",
      "n=60 thresh=3.33e-01 6.28e-02s 3090 + 1642 = 4732 == 4732\n",
      "n=60 thresh=4.44e-01 5.09e-02s 2244 + 1752 = 3996 == 3996\n",
      "n=60 thresh=5.56e-01 3.00e-02s 1450 + 1760 = 3210 == 3210\n",
      "n=60 thresh=6.67e-01 1.31e-02s 846 + 1647 = 2493 == 2493\n",
      "n=60 thresh=7.78e-01 7.54e-03s 358 + 1219 = 1577 == 1577\n",
      "n=60 thresh=8.89e-01 1.82e-03s 94 + 758 = 852 == 852\n",
      "n=60 thresh=1.00e+00 4.67e-04s 0 + 0 = 0 == 0\n",
      "n=100 thresh=0.00e+00 7.12e-01s 19948 + 52 = 20000 == 20000\n",
      "n=100 thresh=1.11e-01 4.73e-01s 15786 + 2024 = 17810 == 17810\n",
      "n=100 thresh=2.22e-01 3.77e-01s 11986 + 3517 = 15503 == 15503\n",
      "n=100 thresh=3.33e-01 2.63e-01s 8924 + 4423 = 13347 == 13347\n",
      "n=100 thresh=4.44e-01 1.73e-01s 6148 + 4942 = 11090 == 11090\n",
      "n=100 thresh=5.56e-01 9.14e-02s 3912 + 4965 = 8877 == 8877\n",
      "n=100 thresh=6.67e-01 4.71e-02s 2192 + 4402 = 6594 == 6594\n",
      "n=100 thresh=7.78e-01 2.28e-02s 1076 + 3295 = 4371 == 4371\n",
      "n=100 thresh=8.89e-01 4.87e-03s 258 + 1956 = 2214 == 2214\n",
      "n=100 thresh=1.00e+00 1.13e-03s 0 + 0 = 0 == 0\n"
     ]
    }
   ],
   "source": [
    "# Unit test for matchEdgeMaps2D.\n",
    "for n in np.round(np.logspace(0,2,10)).astype(int):\n",
    "    for thresh in np.linspace(0,1,10):\n",
    "        print 'n={:d} thresh={:1.2e}'.format(n,thresh),\n",
    "        bmap1 = (np.random.rand(n,n)>thresh).astype(float)\n",
    "        bmap2 = (np.random.rand(n,n)>thresh).astype(float)\n",
    "        t = time.time()\n",
    "        cost,match1,match2 = matchEdgeMaps2D(bmap1,bmap2,.005,100)\n",
    "        print '{:1.2e}s'.format(time.time() - t),\n",
    "        npix = len(bmap1.nonzero()[0]) + len(bmap2.nonzero()[0])\n",
    "        ncor = len(match1.nonzero()[0]) + len(match1.nonzero()[0])\n",
    "        nout = len(((bmap1>0)&(match1==0)).nonzero()[0]) + len(((bmap2>0)&(match2==0)).nonzero()[0])\n",
    "        print ncor,'+',nout,'=',ncor+nout,'==',npix\n",
    "        if ncor+nout!=npix: raise ValueError(\"bug\")\n",
    "        \n",
    "        match1f = match1.ravel()\n",
    "        idx1 = match1f.astype(int)\n",
    "        a1 = match1f.nonzero()[0]\n",
    "        a2 = match1f[idx1.nonzero()[0]].ravel()\n",
    "        a = vstack((a1,a2)).T\n",
    "\n",
    "        match2f = match2.ravel()\n",
    "        idx2 = match2f.astype(int)\n",
    "        b1 = match2f.nonzero()[0]\n",
    "        b2 = match2f[idx2.nonzero()[0]].ravel()\n",
    "        b = vstack((b2,b1)).T\n",
    "        b = b[b[:,0].argsort(),]\n",
    "        if np.sum(a-b): raise ValueError(\"bug\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "n=1 thresh=5.00e-01 0.0001s 2 + 0 = 2 == 2\n",
      "n=1 thresh=6.33e-01 0.0000s 0 + 1 = 1 == 1\n",
      "n=1 thresh=7.67e-01 0.0000s 0 + 0 = 0 == 0\n",
      "n=1 thresh=9.00e-01 0.0000s 0 + 0 = 0 == 0\n",
      "n=4 thresh=5.00e-01 0.0003s 34 + 28 = 62 == 62\n",
      "n=4 thresh=6.33e-01 0.0001s 12 + 38 = 50 == 50\n",
      "n=4 thresh=7.67e-01 0.0001s 6 + 23 = 29 == 29\n",
      "n=4 thresh=9.00e-01 0.0000s 0 + 7 = 7 == 7\n",
      "n=16 thresh=5.00e-01 0.0406s 2018 + 2061 = 4079 == 4079\n",
      "n=16 thresh=6.33e-01 0.0236s 1190 + 1884 = 3074 == 3074\n",
      "n=16 thresh=7.67e-01 0.0070s 434 + 1445 = 1879 == 1879\n",
      "n=16 thresh=9.00e-01 0.0015s 98 + 707 = 805 == 805\n",
      "n=64 thresh=5.00e-01 9.7403s 130884 + 131329 = 262213 == 262213\n",
      "n=64 thresh=6.33e-01 3.4553s 70100 + 121800 = 191900 == 191900\n",
      "n=64 thresh=7.67e-01 1.2439s 28474 + 93913 = 122387 == 122387\n",
      "n=64 thresh=9.00e-01 0.1635s 5222 + 47354 = 52576 == 52576\n"
     ]
    }
   ],
   "source": [
    "# Unit test for correspondVoxels\n",
    "for n in np.round(np.logspace(0,6,4,base=2)).astype(int):\n",
    "    for thresh in np.linspace(.5,.9,4):\n",
    "        print 'n={:d} thresh={:1.2e}'.format(n,thresh),\n",
    "        bmap1 = (np.random.rand(n,n,n)>thresh).astype(float)\n",
    "        bmap2 = (np.random.rand(n,n,n)>thresh).astype(float)\n",
    "\n",
    "        t = time.time()\n",
    "        cost,oc,match1,match2 = correspondVoxels(bmap1,bmap2,.005)\n",
    "        print '{:1.4f}s'.format(time.time() - t),\n",
    "        npix = len(bmap1.nonzero()[0]) + len(bmap2.nonzero()[0])\n",
    "        ncor = len(match1.nonzero()[0]) + len(match1.nonzero()[0])\n",
    "        nout = len(((bmap1>0)&(match1==0)).nonzero()[0]) + len(((bmap2>0)&(match2==0)).nonzero()[0])\n",
    "        print ncor,'+',nout,'=',ncor+nout,'==',npix\n",
    "        if ncor+nout!=npix: raise ValueError(\"bug\")\n",
    "        \n",
    "        match1f = match1.ravel()\n",
    "        idx1 = match1f.astype(int)\n",
    "        a1 = match1f.nonzero()[0]\n",
    "        a2 = match1f[idx1.nonzero()[0]].ravel()\n",
    "        a = vstack((a1,a2)).T\n",
    "\n",
    "        match2f = match2.ravel()\n",
    "        idx2 = match2f.astype(int)\n",
    "        b1 = match2f.nonzero()[0]\n",
    "        b2 = match2f[idx2.nonzero()[0]].ravel()\n",
    "        b = vstack((b2,b1)).T\n",
    "        b = b[b[:,0].argsort(),]\n",
    "        if np.sum(a-b): raise ValueError(\"bug\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "n=1 thresh=5.00e-01 0.0043s 2 + 0 = 2 == 2\n",
      "n=1 thresh=6.33e-01 0.0000s 0 + 0 = 0 == 0\n",
      "n=1 thresh=7.67e-01 0.0000s 0 + 0 = 0 == 0\n",
      "n=1 thresh=9.00e-01 0.0000s 0 + 0 = 0 == 0\n",
      "n=4 thresh=5.00e-01 0.0003s 28 + 37 = 65 == 65\n",
      "n=4 thresh=6.33e-01 0.0002s 14 + 31 = 45 == 45\n",
      "n=4 thresh=7.67e-01 0.0001s 6 + 20 = 26 == 26\n",
      "n=4 thresh=9.00e-01 0.0001s 8 + 15 = 23 == 23\n",
      "n=16 thresh=5.00e-01 0.0417s 2046 + 2050 = 4096 == 4096\n",
      "n=16 thresh=6.33e-01 0.0231s 1130 + 1877 = 3007 == 3007\n",
      "n=16 thresh=7.67e-01 0.0066s 454 + 1484 = 1938 == 1938\n",
      "n=16 thresh=9.00e-01 0.0015s 82 + 730 = 812 == 812\n",
      "n=64 thresh=5.00e-01 31.7409s 130928 + 131151 = 262079 == 262079\n",
      "n=64 thresh=6.33e-01 25.7225s 70730 + 121897 = 192627 == 192627\n",
      "n=64 thresh=7.67e-01 23.0356s 28522 + 93828 = 122350 == 122350\n",
      "n=64 thresh=9.00e-01 9.9525s 5338 + 47387 = 52725 == 52725\n"
     ]
    }
   ],
   "source": [
    "# Unit test for matchEdgeMaps3D\n",
    "for n in np.round(np.logspace(0,6,4,base=2)).astype(int):\n",
    "    for thresh in np.linspace(.5,.9,4):\n",
    "        print 'n={:d} thresh={:1.2e}'.format(n,thresh),\n",
    "        bmap1 = (np.random.rand(n,n,n)>thresh).astype(float)\n",
    "        bmap2 = (np.random.rand(n,n,n)>thresh).astype(float)\n",
    "\n",
    "        t = time.time()\n",
    "        cost,match1,match2 = matchEdgeMaps3D(bmap1,bmap2,.005*np.sqrt(n**3),100)\n",
    "        print '{:1.4f}s'.format(time.time() - t),\n",
    "        npix = len(bmap1.nonzero()[0]) + len(bmap2.nonzero()[0])\n",
    "        ncor = len(match1.nonzero()[0]) + len(match1.nonzero()[0])\n",
    "        nout = len(((bmap1>0)&(match1==0)).nonzero()[0]) + len(((bmap2>0)&(match2==0)).nonzero()[0])\n",
    "        print ncor,'+',nout,'=',ncor+nout,'==',npix\n",
    "        if ncor+nout!=npix: raise ValueError(\"bug\")\n",
    "        \n",
    "        match1f = match1.ravel()\n",
    "        idx1 = match1f.astype(int)\n",
    "        a1 = match1f.nonzero()[0]\n",
    "        a2 = match1f[idx1.nonzero()[0]].ravel()\n",
    "        a = vstack((a1,a2)).T\n",
    "\n",
    "        match2f = match2.ravel()\n",
    "        idx2 = match2f.astype(int)\n",
    "        b1 = match2f.nonzero()[0]\n",
    "        b2 = match2f[idx2.nonzero()[0]].ravel()\n",
    "        b = vstack((b2,b1)).T\n",
    "        b = b[b[:,0].argsort(),]\n",
    "        if np.sum(a-b): raise ValueError(\"bug\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "n=10 thresh=5.00e-01 6.77e-04s 98 48 50\n"
     ]
    }
   ],
   "source": [
    ""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    ""
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2.0
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}